!================Programa que calcula algumas integrais=================

PROGRAM exerB

!=======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis=========================

      INTEGER :: num_N, I, J, min_err
      INTEGER, ALLOCATABLE, DIMENSION(:) :: N_values
      REAL(8) :: f, prim_f, prim_value, mid
      REAL(8), ALLOCATABLE, DIMENSION(:) :: h
      REAL(8), ALLOCATABLE, DIMENSION(:,:) :: tab
      !Limites de integração (OBS: a deve ser o menor limite)
      REAL(8), PARAMETER :: a = 0.0d0, b = 1.0d0
      !Formato a ser usado no write que escreve no arquivo
      CHARACTER(LEN=*), PARAMETER :: FMT1 = "(5ES16.7)"

!==========================Leitura de dados=============================

      OPEN(UNIT=10, FILE='tab2_in.dat', STATUS='old')
      OPEN(UNIT=20, FILE='tab2_out.dat', STATUS='unknown')

      !Lê o número de valores de N
      READ(10,*) num_N

      !Aloca os vetores que contém os valores de N e h
      ALLOCATE(N_values(num_N),h(num_N),tab(5,num_N))

      !Lê os valores de N
      READ(10,*) N_values

      tab = 0.0d0

!=========================Gera os valores de h==========================

      DO I=1,num_N
            h(I) = (b - a)/N_values(I)
      END DO

!====Escreve N e h na tabela e também calcula por Trapézio e Simpson====

      !Calcula o valor da integral pela função primitiva
      prim_value = prim_f(b) - prim_f(a)

      DO I=1,num_N
            !Escreve N na tabela
            tab(1,I) = N_values(I)
            !Escreve h na tabela
            tab(2,I) = h(I)
            !Inicializa as variáveis a serem usadas no DO WHILE
            J = 1
            mid = 0.0d0
            tab(3,I) = 0.0d0
            tab(4,I) = 0.0d0
            !Calculo das integrais
            DO WHILE (mid <= (b-h(I)))
                  !Calcula o meio do intervalo a ser usado como
                  !como referencia no cálculo da área
                  mid = a + J*h(I)
                  !Regra do Trapézio
                  tab(3,I) = tab(3,I) + (h(I)/2.0d0)*(f(mid+h(I))+&
                  &2.0d0*f(mid)+f(mid-h(I)))
                  !Regra de Simpson
                  tab(4,I) = tab(4,I) + (h(I)/3.0d0)*(f(mid+h(I))+&
                  &4.0d0*f(mid)+f(mid-h(I)))
                  J = J + 2
            END DO
            tab(3,I) = ABS(tab(3,I) - prim_value)
            tab(4,I) = ABS(tab(4,I) - prim_value)
      END DO

!==============================Regra de Bode============================

      DO I=1,num_N
            !Inicializa variáveis a serem usadas no DO WHILE
            J = 0
            mid = 0.0d0 
            tab(5,I) = 0.0d0
            DO WHILE (mid <= (b-4*h(I)))
                  mid = a + J*4*h(I)
                  !Regra de Bode
                  tab(5,I) = tab(5,I) + (2.0d0*h(I)/45.0d0)*(7.0d0*&
                  &f(mid)+32.0d0*f(mid+h(I))+12.0d0*f(mid+2.0d0*h(I))&
                  &+32.0d0*f(mid+3.0d0*h(I))+7.0d0*f(mid+4.0d0*h(I)))
                  J = J + 1
            END DO
            tab(5,I) = ABS(tab(5,I) - prim_value)
      END DO

!=========================Escreve no arquivo============================

      WRITE(20,FMT1) tab

!==========================Acha o melhor N==============================

      DO I=3,5
            DO J=1,(num_N-1)
                  IF (tab(I,J+1) < tab(I,J)) THEN
                        min_err = J+1
                  END IF
            END DO
            WRITE(*,*) "No", I-2, "º caso, o melhor N é:", &
            &tab(1,min_err)
      END DO

!==========================Fecha os arquivos============================

      CLOSE(UNIT=10)
      CLOSE(UNIT=20)

!=======================================================================

END PROGRAM exerB

!=======================================================================


!========================Subroutines e funções==========================

!Calcula o valor de f no ponto x
REAL(8) FUNCTION f(x)
      IMPLICIT NONE

      REAL(8) :: x

      f = DEXP(x/4.0d0)*DSIN(x)
END FUNCTION f

!Calcula o valor da primitiva de f no ponto x
REAL(8) FUNCTION prim_f(x)
      IMPLICIT NONE

      REAL(8) :: x

      prim_f = (-4.0d0/17.0d0)*DEXP(x/4.0d0)*(4.0d0*DCOS(x)-DSIN(x))
END FUNCTION prim_f
